Billiards and Tilts Analysis Notebook

library(plotly)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

0.1 Load all of the throw data

# Note: this notebook right now REQUIRES that the analysis notebook be run to
# make the fits, or that the data folder contains the fits

# load omnibus dataframe
omnibus_df <- read_delim("../data/processed/omnibus/omnibus_raw.csv",
  delim = ",",
  col_types = cols(
    .default = col_double(),
    type = col_factor(),
    ppid = col_factor(),
    exp_label = col_factor(),
    experiment = col_factor(),
    hand = col_factor(),
    camera_tilt = col_factor(),
    surface_tilt = col_factor(),
    target = col_factor(),
    test_type = col_factor(),
    phase = col_factor(),
    learning_and_decay_curves = col_factor(),
    prior_anim = col_factor(),
    baseline_block = col_logical(),
    alt_washout_block = col_logical(),
    alt_all_washout_block = col_logical(),
    task_type = col_factor(),
    surface = col_factor(),
    anim_type = col_factor()
  )
) %>% # filter out practice blocks
  filter(
    block_num > 4,
    !(ppid %in% c(2, 81, 82, 83))
  )

# load exponential decay fits
init_learning_rates <- read_csv(
  "../data/processed/learning_rate_df.csv",
  col_types = cols(
    .default = col_double(),
    experiment = col_factor(),
    phase = col_factor()
  )
) %>%
  mutate(experiment = factor(
    experiment,
    levels = c(
      "accel_cued_tilt", "accel_uncued",
      "curved_cued_tilt", "curved_uncued",
      "rot30_cued_tilt", "rot30_uncued",
      "rot15_cued_tilt", "rot15_uncued",
      "a_ball_roll_animate_surface"
    )
  ))

init_learning_rates_3par <- read_csv(
  "../data/processed/learning_rate_df_3par.csv",
  col_types = cols(
    .default = col_double(),
    experiment = col_factor(),
    phase = col_factor()
  )
) %>%
  mutate(experiment = factor(
    experiment,
    levels = c(
      "accel_cued_tilt", "accel_uncued",
      "curved_cued_tilt", "curved_uncued",
      "rot30_cued_tilt", "rot30_uncued",
      "rot15_cued_tilt", "rot15_uncued",
      "a_ball_roll_animate_surface"
    )
  ))

anim_learning_rates <- read_csv(
  "../data/processed/learning_rate_df_anim.csv",
  col_types = cols(
    .default = col_double(),
    ppid = col_factor(),
    experiment = col_factor(),
    phase = col_factor()
  )
) %>%
  mutate(experiment = factor(
    experiment,
    levels = c(
      "rot30_cued_tilt", "rot30_uncued",
      "half_anim", "full_anim"
    )
  ))

alt_washout_rates <- read_csv(
  "../data/processed/exp_fits_alt_washout_curves.csv",
  col_types = cols(
    .default = col_double(),
    experiment = col_factor(),
    phase = col_factor()
    )
  ) %>%
  mutate(experiment = factor(
    experiment, 
    levels = c(
      "accel_cued_tilt", "accel_uncued",
      "curved_cued_tilt", "curved_uncued",
      "rot30_cued_tilt", "rot30_uncued",
      "rot15_cued_tilt", "rot15_uncued",
      "a_ball_roll_animate_surface"
      )
    )
  )

fitted_learning_rate_df <- read_csv(
  "../data/processed/fit_curves/fitted_learning_rate_df.csv",
  col_types = cols(
    .default = col_double(),
    experiment = col_factor(),
    phase = col_factor()
  )
) %>%
  mutate(experiment = factor(
    experiment,
    levels = c(
      "accel_cued_tilt", "accel_uncued",
      "curved_cued_tilt", "curved_uncued",
      "rot30_cued_tilt", "rot30_uncued",
      "rot15_cued_tilt", "rot15_uncued",
      "a_ball_roll_animate_surface"
    )
  ))

fitted_alt_washout_df <- read_csv(
  "../data/processed/fit_curves/fitted_exp_fits_alt_washout_curves.csv",
  col_types = cols(
    .default = col_double(),
    experiment = col_factor(),
    phase = col_factor()
  )
) %>%
  mutate(experiment = factor(
    experiment,
    levels = c(
      "accel_cued_tilt", "accel_uncued",
      "curved_cued_tilt", "curved_uncued",
      "rot30_cued_tilt", "rot30_uncued",
      "rot15_cued_tilt", "rot15_uncued",
      "a_ball_roll_animate_surface"
    )
  ))

# use flick_direction_x, y, and z to get magnitude of flick direction
omnibus_df$magnitude <- sqrt(omnibus_df$flick_direction_x^2 + omnibus_df$flick_direction_y^2 + omnibus_df$flick_direction_z^2)

1 Success Manifolds - Accelerated Groups

For these plots, we will use data only from the training phase (transfer too?). Baseline phases and washout phases are the same for all groups.

data_per_trial <- omnibus_df %>%
  filter(
    phase == "training",
    experiment %in% c("accel_uncued", "accel_cued_tilt")
  )

# Note for figure: facet by target

# set up the plot
p <- data_per_trial %>%
  ggplot(aes(
    x = raw_throw_deviation,
    y = magnitude,
    color = raw_error_size
  )) +
  # vertical line at 0
  # geom_vline(xintercept = 0, linetype = "dashed", color = "#CCCCCC") +
  # colour palette
  scale_colour_gradient(
    low = "#7eaee2", high = "#001a35",
    limits = c(.05, .70), na.value = "white" # Note: max error value is ~68, target radius is 5
  ) +
  # add data points
  geom_point(size = 1.3, alpha = 0.7, stroke = 0) +
  # facet
  facet_wrap(~ fct_relevel(target, "96", "92", "88", "84"), nrow = 1) +
  # theme
  theme_classic() +
  labs(
    x = NULL,
    y = "Throw Speed (m/s)"
  ) +
  theme(
    panel.background = element_rect(
      fill = "#000000",
      colour = "#000000"
    ),
    legend.position = "none",
    text = element_text(size = text_size),
    strip.background = element_blank(), # facet labels
    strip.text.x = element_blank()
  ) +
  scale_y_continuous(limits = c(1.7, 2.2)) +
  scale_x_continuous(limits = c(-60, 10))

# save
if (save_plots) {
  ggsave(
    p,
    filename = "../data/figs/paper_figs/fig3_manifolds_accel.pdf", device = "pdf",
    height = 1.7, width = 7.2
  )
}

ggplotly(p)

# p

data_per_trial <- omnibus_df %>%
  filter(
    phase == "training",
    experiment %in% c("rot30_uncued", "rot30_cued_tilt")
  )

# Note for figure: facet by target

# set up the plot
p <- data_per_trial %>%
  ggplot(aes(
    x = raw_throw_deviation,
    y = magnitude,
    color = raw_error_size
    # shape = ppid
  )) +
  # vertical line at 0
  geom_vline(xintercept = -30, linetype = "dashed", color = "#CCCCCC") +
  # colour palette
  scale_colour_gradient(
    low = "#7eaee2", high = "#001a35",
    limits = c(.05, .70), na.value = "white" # Note: max error value is ~68, target radius is 5
  ) +
  # add data points
  geom_point(size = 1.3, alpha = 0.7, stroke = 0) +
  # facet
  facet_wrap(~ fct_relevel(target, "96", "92", "88", "84"), nrow = 1) +
  # theme
  theme_classic() +
  labs(
    x = NULL,
    y = "Throw Speed (m/s)"
  ) +
  theme(
    panel.background = element_rect(
      fill = "#000000",
      colour = "#000000"
    ),
    legend.position = "none",
    text = element_text(size = text_size),
    strip.background = element_blank(), # facet labels
    strip.text.x = element_blank()
  ) +
  scale_y_continuous(limits = c(1.7, 2.2)) +
  scale_x_continuous(limits = c(-60, 10))

# save
if (save_plots) {
  ggsave(
    p,
    filename = "../data/figs/paper_figs/fig3_manifolds_rotated.pdf", device = "pdf",
    height = 1.7, width = 7.2
  )
}

# ggplotly(p)
p

data_per_trial <- omnibus_df %>%
  filter(
    phase == "training",
    experiment %in% c("curved_cued_tilt", "curved_uncued")
  )

# Note for figure: facet by target

# set up the plot
p <- data_per_trial %>%
  ggplot(aes(
    x = raw_throw_deviation,
    y = magnitude,
    color = raw_error_size
    # shape = ppid
  )) +
  # vertical line at 30
  geom_vline(xintercept = -30, linetype = "dashed", color = "#CCCCCC") +
  # colour palette
  scale_colour_gradient(
    low = "#7eaee2", high = "#001a35",
    limits = c(.05, .70), na.value = "white" # Note: max error value is ~68, target radius is 5
  ) +
  # add data points
  geom_point(size = 1.3, alpha = 0.7, stroke = 0) +
  # facet
  facet_wrap(~ fct_relevel(target, "96", "92", "88", "84"), nrow = 1) +
  # theme
  theme_classic() +
  labs(
    x = NULL,
    y = "Throw Speed (m/s)"
  ) +
  theme(
    panel.background = element_rect(
      fill = "#000000",
      colour = "#000000"
    ),
    legend.position = "none",
    text = element_text(size = text_size),
    strip.background = element_blank(), # facet labels
    strip.text.x = element_blank()
  ) +
  scale_y_continuous(limits = c(1.7, 2.2)) +
  scale_x_continuous(limits = c(-60, 10))

# save
if (save_plots) {
  ggsave(
    p,
    filename = "../data/figs/paper_figs/fig3_manifolds_curved.pdf", device = "pdf",
    height = 1.7, width = 7.2
  )
}

# ggplotly(p)
p

num_exps <- 4

1.1 Normalized Training Phases

# rest of the exps
data_per_group <- omnibus_df %>%
  filter(
    experiment %in% c(
      "rot30_cued_tilt", "rot30_uncued",
      "curved_cued_tilt", "curved_uncued"
    ),
    trial_num >= 214,
    phase == "training",
    learning_and_decay_curves == 1
  ) %>%
  group_by(experiment, phase, test_type, trial_num) %>%
  summarise(
    mean_deviation = mean(norm_throw_deviation),
    ci_deviation = vector_confint(norm_throw_deviation),
    .groups = "drop"
  ) %>%
  arrange(experiment, trial_num)

# add a dummy column with repeating sequence
# NOTE: this can't be combined with above since we are using nrow
data_per_group <- data_per_group %>%
  mutate(dummy_x = rep(1:(nrow(data_per_group) / num_exps),
    length.out = nrow(data_per_group)
  ))

# set up plot
p <- data_per_group %>%
  ggplot(
    aes(
      x = dummy_x, y = mean_deviation, colour = experiment
    )
  ) +
  theme_classic() +
  theme(
    legend.position = "none",
    text = element_text(size = text_size)
  ) +
  labs(
    x = "Trial",
    y = "Normalized Throw Angle"
  ) +
  # add horizontal lines
  geom_hline(
    yintercept = c(0, 1), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "solid"
  ) +
  geom_hline(
    yintercept = c(0.5, 1.5), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "dashed"
  ) +
  # set colour palette
  scale_colour_manual(values = pallete_list) +
  scale_fill_manual(values = pallete_list) + # scale axes
  scale_y_continuous(
    breaks = c(0, 0.5, 1, 1.5),
    labels = c(0, 0.5, 1, 1.5),
    limits = c(-0.2, 1.5)
  ) +
  scale_x_continuous(
    breaks = curve_x_breaks,
    labels = curve_x_breaks
  )

# add confidence intervals and data points
for (uniq_phase in unique(data_per_group$phase)) {
  # get the data for this block
  to_plot_data <- filter(data_per_group, phase == uniq_phase)

  p <- p + geom_ribbon(
    data = to_plot_data,
    aes(
      ymin = mean_deviation - ci_deviation,
      ymax = mean_deviation + ci_deviation,
      fill = experiment
    ), colour = NA, alpha = 0.3
  ) + geom_line(
    data = to_plot_data
  )
}

# save
if (save_plots) {
  ggsave(
    p,
    filename = "../data/figs/paper_figs/fig3_exps_normalized_training.pdf", device = "pdf",
    height = 2, width = curve_fig_width
  )
}

p

fitted_data <- fitted_learning_rate_df %>%
  filter(
    experiment %in% c(
      "rot30_cued_tilt", "rot30_uncued",
      "curved_cued_tilt", "curved_uncued"
    ),
    phase == "training"
  ) %>%
  mutate(
    pert = factor(str_split_fixed(experiment, "_", 2)[, 1]),
    cue = factor(str_split_fixed(experiment, "_", 2)[, 2])
  )

p <- fitted_data %>%
  ggplot(
    aes(x = x + 1, y = y, colour = experiment)
  ) +
  theme_classic() +
  theme(
    legend.position = "none",
    text = element_text(size = text_size)
  ) +
  labs(
    x = "Trial",
    y = "Normalized Throw Angle"
  ) +
  # add horizontal lines
  geom_hline(
    yintercept = c(0, 1), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "solid"
  ) +
  geom_hline(
    yintercept = c(0.5, 1.5), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "dashed"
  ) +
  # set colour palette
  scale_colour_manual(values = pallete_list) +
  scale_fill_manual(values = pallete_list) + # scale axes
  scale_y_continuous(
    breaks = c(0, 0.5, 1, 1.5),
    labels = c(0, 0.5, 1, 1.5),
    limits = c(-0.2, 1.5)
  ) +
  scale_x_continuous(
    breaks = curve_x_breaks,
    labels = curve_x_breaks
  ) +
  geom_line()

# save
if (save_plots) {
  ggsave(
    p,
    filename = "../data/figs/paper_figs/fig3_curved_training_fitted_only.pdf", device = "pdf",
    height = 2, width = curve_fig_width
  )
}

p

1.2 Training - Fit Learning Rates

data_per_ppt <- init_learning_rates %>%
  filter(
    experiment %in% c(
      "rot30_cued_tilt", "rot30_uncued",
      "curved_cued_tilt", "curved_uncued"
    ),
    phase == "training"
  ) %>%
  mutate(
    pert = factor(str_split_fixed(experiment, "_", 2)[, 1]),
    cue = factor(str_split_fixed(experiment, "_", 2)[, 2]),
    ppid = factor(ppid)
  )

data_per_group <- data_per_ppt %>%
  group_by(experiment, phase) %>%
  summarise(
    mean_learning_rate = mean(exp_fit_lambda),
    ci_learning_rate = vector_confint(exp_fit_lambda),
    mean_high = mean(exp_fit_N0),
    ci_high = vector_confint(exp_fit_N0),
    n = n(),
    .groups = "drop"
  )

p <- data_per_group %>%
  ggplot(
    aes(x = experiment, y = mean_learning_rate, colour = experiment)
  ) +
  theme_classic() +
  labs(
    x = per_group_x_labs,
    y = "Learning Rate"
  ) +
  # remove all x axis labels
  theme(
    legend.position = "none",
    text = element_text(size = text_size)
  ) +
  # colour legend settings
  guides(
    colour = guide_legend(override.aes = list(alpha = 1))
  ) +
  # set colour palette
  scale_colour_manual(values = pallete_list) +
  scale_fill_manual(
    values = pallete_list
  ) +
  # add horizontal lines
  geom_hline(
    yintercept = c(0, 1), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "solid"
  ) +
  geom_hline(
    yintercept = c(0.25, 0.5, 0.75), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "dashed"
  ) +
  # y-axis scale
  scale_y_continuous(
    # limits = c(-35, 5),
    breaks = c(0, 0.25, 0.5, 0.75, 1),
    labels = c(0, 0.25, 0.5, 0.75, 1)
  ) +
  scale_x_discrete(
    labels = c(0, 0, 0, 0)
  ) +
  # add data points
  geom_beeswarm(
    data = data_per_ppt,
    aes(
      y = exp_fit_lambda
    ),
    alpha = 0.1,
    size = 1
  ) +
  geom_linerange(aes(
    ymin = mean_learning_rate - ci_learning_rate,
    ymax = mean_learning_rate + ci_learning_rate
  ), alpha = 0.5, lwd = 2) +
  geom_point()

# save
if (save_plots) {
  ggsave(
    p,
    filename = "../data/figs/paper_figs/fig3_training_learning_rates.pdf", device = "pdf",
    height = 2, width = per_group_fig_width
  )
}

p

1.2.1 ANOVA - 3 levels

(temp_ANOVA <- aov_car(
  exp_fit_lambda ~ pert * cue + Error(ppid),
  data_per_ppt,
  include_aov = TRUE
))

Contrasts set to contr.sum for the following variables: pert, cue

Anova Table (Type 3 tests)

Response: exp_fit_lambda
    Effect    df  MSE         F  ges p.value
1     pert 1, 77 0.06 12.88 *** .143   <.001
2      cue 1, 77 0.06   7.60 ** .090    .007
3 pert:cue 1, 77 0.06    3.11 + .039    .082
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

print("Bayes anovas and bayes factors:")

[1] "Bayes anovas and bayes factors:"

print("")

[1] ""

bf <- anovaBF(exp_fit_lambda ~ pert * cue, data = data.frame(data_per_ppt), progress = FALSE)
print(bf)

Bayes factor analysis
--------------
[1] pert                  : 33.67806 ±0%
[2] cue                   : 4.088725 ±0.01%
[3] pert + cue            : 190.7611 ±0.64%
[4] pert + cue + pert:cue : 201.0003 ±0.86%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

print(bayesfactor_inclusion(bf))

Inclusion Bayes Factors (Model Averaged)

         P(prior) P(posterior) Inclusion BF
pert         0.60         0.99        55.74
cue          0.60         0.92         7.61
cue:pert     0.20         0.47         3.50

* Compared among: all models
*    Priors odds: uniform-equal

# Note: default multiplicity correction for pairs summary is Tukey HSD
print("----------------------------- EMMs -----------------------------")

[1] "----------------------------- EMMs -----------------------------"

emms_obj <- emmeans(temp_ANOVA, "pert")

NOTE: Results may be misleading due to involvement in interactions

print(pwpm(emms_obj, diffs = FALSE))

        curved   rot30
curved [0.527]  0.0006
rot30          [0.329]

Row and column labels: pert
Upper triangle: P values 
Diagonal: [Estimates] (emmean) 

print("----------------------------- EMMs -----------------------------")

[1] "----------------------------- EMMs -----------------------------"

emms_obj <- emmeans(temp_ANOVA, "cue")

NOTE: Results may be misleading due to involvement in interactions

print(pwpm(emms_obj, diffs = FALSE))

          cued_tilt  uncued
cued_tilt   [0.504]  0.0073
uncued              [0.352]

Row and column labels: cue
Upper triangle: P values 
Diagonal: [Estimates] (emmean) 

# Effect size is Cohen's D
eff_size(emms_obj, edf = temp_ANOVA$anova_table$`den Df`[1], sigma = sigma(temp_ANOVA$lm))

 contrast           effect.size    SE df lower.CL upper.CL
 cued_tilt - uncued       0.613 0.228 77    0.159     1.07

Results are averaged over the levels of: pert 
sigma used for effect sizes: 0.2484 
Confidence level used: 0.95 

print("----------------------------- Full Tukey HSD -----------------------------")

[1] "----------------------------- Full Tukey HSD -----------------------------"

print(TukeyHSD(temp_ANOVA$aov))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = tmp_formula, data = dat.ret, contrasts = contrasts)

$pert
                   diff        lwr         upr     p adj
rot30-curved -0.2012278 -0.3111615 -0.09129398 0.0004837

$cue
                      diff        lwr         upr     p adj
uncued-cued_tilt -0.153335 -0.2632687 -0.04340118 0.0068789

$`pert:cue`
                                        diff        lwr         upr     p adj
rot30:cued_tilt-curved:cued_tilt -0.10082813 -0.3071210  0.10546473 0.5761442
curved:uncued-curved:cued_tilt   -0.05482835 -0.2611212  0.15146451 0.8975670
rot30:uncued-curved:cued_tilt    -0.35036936 -0.5541916 -0.14654716 0.0001305
curved:uncued-rot30:cued_tilt     0.04599978 -0.1602931  0.25229263 0.9361528
rot30:uncued-rot30:cued_tilt     -0.24954123 -0.4533634 -0.04571903 0.0100884
rot30:uncued-curved:uncued       -0.29554100 -0.4993632 -0.09171880 0.0015670

print(eff_size(emmeans(temp_ANOVA$aov, c("pert", "cue")),
  edf = temp_ANOVA$anova_table$`den Df`[1],
  sigma = sigma(temp_ANOVA$lm)
))

 contrast                           effect.size    SE df lower.CL upper.CL
 curved cued_tilt - rot30 cued_tilt       0.406 0.318 77   -0.227    1.039
 curved cued_tilt - curved uncued         0.221 0.317 77   -0.410    0.851
 curved cued_tilt - rot30 uncued          1.410 0.332 77    0.748    2.072
 rot30 cued_tilt - curved uncued         -0.185 0.317 77   -0.816    0.445
 rot30 cued_tilt - rot30 uncued           1.005 0.323 77    0.362    1.647
 curved uncued - rot30 uncued             1.190 0.327 77    0.539    1.840

sigma used for effect sizes: 0.2484 
Confidence level used: 0.95 

1.3 Training - Fit High Points (Asymptotes)

p <- data_per_group %>%
  ggplot(
    aes(x = experiment, y = mean_high, colour = experiment)
  ) +
  theme_classic() +
  labs(
    x = per_group_x_labs,
    y = "Asymptote of Adaptation"
  ) +
  # for the colour legend, only show the first 7
  # Note this doesn't work for the plotly plot
  guides(
    colour = guide_legend(override.aes = list(alpha = 1))
  ) +
  # set colour palette
  scale_colour_manual(values = pallete_list) +
  scale_fill_manual(values = pallete_list) +
  # add horizontal lines
  geom_hline(
    yintercept = c(1), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "solid"
  ) +
  geom_hline(
    yintercept = c(0.5, 1.5), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "dashed"
  ) +
  # y-axis scale
  scale_y_continuous(
    limits = c(0.5, 1.75),
    breaks = c(0.5, 0.75, 1, 1.25, 1.5),
    labels = c(0.5, 0.75, 1, 1.25, 1.5)
  ) +
  scale_x_discrete(
    labels = c(0, 0, 0, 0)
  ) +
  # remove all x axis labels
  theme(
    legend.position = "none",
    text = element_text(size = text_size)
  ) +
  # add data points
  geom_beeswarm(
    data = data_per_ppt,
    aes(
      y = exp_fit_N0
    ),
    alpha = 0.1,
    size = 1
  ) +
  geom_linerange(aes(
    ymin = mean_high - ci_high,
    ymax = mean_high + ci_high
  ), alpha = 0.5, lwd = 2) +
  geom_point()

# save
if (save_plots) {
  ggsave(
    p,
    filename = "../data/figs/paper_figs/fig3_training_asymptotes.pdf", device = "pdf",
    height = 2, width = per_group_fig_width
  )
}

p

1.3.1 ANOVA - 3 levels

(temp_ANOVA <- aov_car(
  exp_fit_N0 ~ pert * cue + Error(ppid),
  data_per_ppt,
  include_aov = TRUE
))

Contrasts set to contr.sum for the following variables: pert, cue

Anova Table (Type 3 tests)

Response: exp_fit_N0
    Effect    df  MSE    F  ges p.value
1     pert 1, 77 0.02 1.57 .020    .214
2      cue 1, 77 0.02 0.25 .003    .621
3 pert:cue 1, 77 0.02 0.31 .004    .582
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

print("Bayes anovas and bayes factors:")

[1] "Bayes anovas and bayes factors:"

print("")

[1] ""

bf <- anovaBF(exp_fit_N0 ~ pert * cue, data = data.frame(data_per_ppt), progress = FALSE)
print(bf)

Bayes factor analysis
--------------
[1] pert                  : 0.4618912  ±0.01%
[2] cue                   : 0.2583483  ±0.02%
[3] pert + cue            : 0.1183708  ±0.99%
[4] pert + cue + pert:cue : 0.04131614 ±1.34%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

print(bayesfactor_inclusion(bf))

Inclusion Bayes Factors (Model Averaged)

         P(prior) P(posterior) Inclusion BF
pert         0.60         0.33        0.329
cue          0.60         0.22        0.191
cue:pert     0.20         0.02        0.090

* Compared among: all models
*    Priors odds: uniform-equal

# Note: default multiplicity correction for pairs summary is Tukey HSD
print("----------------------------- EMMs -----------------------------")

[1] "----------------------------- EMMs -----------------------------"

emms_obj <- emmeans(temp_ANOVA, "pert")

NOTE: Results may be misleading due to involvement in interactions

print(pwpm(emms_obj, diffs = FALSE))

       curved  rot30
curved [1.08] 0.2142
rot30         [1.03]

Row and column labels: pert
Upper triangle: P values 
Diagonal: [Estimates] (emmean) 

# Effect size is Cohen's D
eff_size(emms_obj, edf = temp_ANOVA$anova_table$`den Df`[1], sigma = sigma(temp_ANOVA$lm))

 contrast       effect.size    SE df lower.CL upper.CL
 curved - rot30       0.278 0.223 77   -0.166    0.723

Results are averaged over the levels of: cue 
sigma used for effect sizes: 0.1479 
Confidence level used: 0.95 

print("----------------------------- Full Tukey HSD -----------------------------")

[1] "----------------------------- Full Tukey HSD -----------------------------"

print(TukeyHSD(temp_ANOVA$aov))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = tmp_formula, data = dat.ret, contrasts = contrasts)

$pert
                    diff        lwr        upr     p adj
rot30-curved -0.04113983 -0.1065787 0.02429902 0.2144126

$cue
                        diff         lwr        upr     p adj
uncued-cued_tilt -0.01610158 -0.08154043 0.04933727 0.6255562

$`pert:cue`
                                         diff         lwr        upr     p adj
rot30:cued_tilt-curved:cued_tilt -0.059316474 -0.18211375 0.06348080 0.5856122
curved:uncued-curved:cued_tilt   -0.034477005 -0.15727428 0.08832027 0.8817483
rot30:uncued-curved:cued_tilt    -0.057484866 -0.17881146 0.06384173 0.6009106
curved:uncued-rot30:cued_tilt     0.024839469 -0.09795781 0.14763675 0.9512295
rot30:uncued-rot30:cued_tilt      0.001831608 -0.11949499 0.12315821 0.9999774
rot30:uncued-curved:uncued       -0.023007861 -0.14433446 0.09831874 0.9592907

print(eff_size(emmeans(temp_ANOVA$aov, c("pert", "cue")),
  edf = temp_ANOVA$anova_table$`den Df`[1],
  sigma = sigma(temp_ANOVA$lm)
))

 contrast                           effect.size    SE df lower.CL upper.CL
 curved cued_tilt - rot30 cued_tilt      0.4011 0.318 77   -0.232    1.034
 curved cued_tilt - curved uncued        0.2332 0.317 77   -0.398    0.864
 curved cued_tilt - rot30 uncued         0.3887 0.314 77   -0.237    1.014
 rot30 cued_tilt - curved uncued        -0.1680 0.317 77   -0.798    0.462
 rot30 cued_tilt - rot30 uncued         -0.0124 0.312 77   -0.635    0.610
 curved uncued - rot30 uncued            0.1556 0.313 77   -0.467    0.778

sigma used for effect sizes: 0.1479 
Confidence level used: 0.95 

1.4 Normalized Washout Phases

# rest of the exps
data_per_group <- omnibus_df %>%
  filter(
    experiment %in% c(
      "rot30_cued_tilt", "rot30_uncued",
      "curved_cued_tilt", "curved_uncued"
    ),
    trial_num >= 214,
    phase == "washout",
    learning_and_decay_curves == 1
  ) %>%
  group_by(experiment, phase, trial_num) %>%
  summarise(
    mean_deviation = mean(norm_throw_deviation),
    ci_deviation = vector_confint(norm_throw_deviation),
    .groups = "drop"
  ) %>%
  arrange(experiment, trial_num)

# add a dummy column with repeating sequence
# NOTE: this can't be combined with above since we are using nrow
data_per_group <- data_per_group %>%
  mutate(dummy_x = rep(1:(nrow(data_per_group) / num_exps),
    length.out = nrow(data_per_group)
  ))

# set up plot
p <- data_per_group %>%
  ggplot(
    aes(
      x = dummy_x, y = mean_deviation, colour = experiment
    )
  ) +
  theme_classic() +
  theme(
    legend.position = "none",
    text = element_text(size = text_size)
  ) +
  labs(
    x = "Trial",
    y = "Normalized Throw Angle"
  ) +
  # add horizontal lines
  geom_hline(
    yintercept = c(0, 1), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "solid"
  ) +
  geom_hline(
    yintercept = c(0.5, 1.5), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "dashed"
  ) +
  # set colour palette
  scale_colour_manual(values = pallete_list) +
  scale_fill_manual(values = pallete_list) + # scale axes
  scale_y_continuous(
    breaks = c(0, 0.5, 1, 1.5),
    labels = c(0, 0.5, 1, 1.5)
  ) +
  scale_x_continuous(
    breaks = curve_x_breaks,
    labels = curve_x_breaks
  )

# add confidence intervals and data points
for (uniq_phase in unique(data_per_group$phase)) {
  # get the data for this block
  to_plot_data <- filter(data_per_group, phase == uniq_phase)

  p <- p + geom_ribbon(
    data = to_plot_data,
    aes(
      ymin = mean_deviation - ci_deviation,
      ymax = mean_deviation + ci_deviation,
      fill = experiment
    ), colour = NA, alpha = 0.3
  ) + geom_line(
    data = to_plot_data
  )
}

# save
if (save_plots) {
  ggsave(
    p,
    filename = "../data/figs/paper_figs/fig3_normalized_washout.pdf", device = "pdf",
    height = 2, width = 2.5
  )
}

p

1.5 Washout - Fit Learning Rates

data_per_ppt <- init_learning_rates %>%
  filter(
    experiment %in% c(
      "rot30_cued_tilt", "rot30_uncued",
      "curved_cued_tilt", "curved_uncued"
    ),
    phase == "washout"
  ) %>%
  mutate(
    pert = factor(str_split_fixed(experiment, "_", 2)[, 1]),
    cue = factor(str_split_fixed(experiment, "_", 2)[, 2]),
    ppid = factor(ppid)
  )

data_per_group <- data_per_ppt %>%
  group_by(experiment, phase) %>%
  summarise(
    mean_learning_rate = mean(exp_fit_lambda),
    ci_learning_rate = vector_confint(exp_fit_lambda),
    mean_high = mean(exp_fit_N0),
    ci_high = vector_confint(exp_fit_N0),
    .groups = "drop"
  )

p <- data_per_group %>%
  ggplot(
    aes(x = experiment, y = mean_learning_rate, colour = experiment)
  ) +
  theme_classic() +
  labs(
    x = per_group_x_labs,
    y = "Decay Rate"
  ) +
  # remove all x axis labels
  theme(
    legend.position = "none",
    text = element_text(size = text_size)
  ) +
  # colour legend settings
  guides(
    colour = guide_legend(override.aes = list(alpha = 1))
  ) +
  # set colour palette
  scale_colour_manual(values = pallete_list) +
  scale_fill_manual(
    values = pallete_list
  ) +
  # add horizontal lines
  geom_hline(
    yintercept = c(0, 1), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "solid"
  ) +
  geom_hline(
    yintercept = c(0.25, 0.5, 0.75), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "dashed"
  ) +
  # y-axis scale
  scale_y_continuous(
    # limits = c(-35, 5),
    breaks = c(0, 0.25, 0.5, 0.75, 1),
    labels = c(0, 0.25, 0.5, 0.75, 1)
  ) +
  scale_x_discrete(
    labels = c(0, 0, 0, 0)
  ) +
  # add data points
  geom_beeswarm(
    data = data_per_ppt,
    aes(
      y = exp_fit_lambda
    ),
    alpha = 0.1,
    size = 1
  ) +
  geom_linerange(aes(
    ymin = mean_learning_rate - ci_learning_rate,
    ymax = mean_learning_rate + ci_learning_rate
  ), alpha = 0.5, lwd = 2) +
  geom_point()

# save
if (save_plots) {
  ggsave(
    p,
    filename = "../data/figs/paper_figs/fig3_washout_learning_rates.pdf", device = "pdf",
    height = 2, width = per_group_fig_width
  )
}

p

1.5.1 ANOVA - 3 levels

(temp_ANOVA <- aov_car(
  exp_fit_lambda ~ pert * cue + Error(ppid),
  data_per_ppt,
  include_aov = TRUE
))

Contrasts set to contr.sum for the following variables: pert, cue

Anova Table (Type 3 tests)

Response: exp_fit_lambda
    Effect    df  MSE         F  ges p.value
1     pert 1, 77 0.07 18.48 *** .194   <.001
2      cue 1, 77 0.07      1.07 .014    .304
3 pert:cue 1, 77 0.07      0.47 .006    .496
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

print("Bayes anovas and bayes factors:")

[1] "Bayes anovas and bayes factors:"

print("")

[1] ""

bf <- anovaBF(exp_fit_lambda ~ pert * cue, data = data.frame(data_per_ppt), progress = FALSE)
print(bf)

Bayes factor analysis
--------------
[1] pert                  : 433.4347  ±0%
[2] cue                   : 0.3244278 ±0.02%
[3] pert + cue            : 150.7467  ±2%
[4] pert + cue + pert:cue : 56.55872  ±1.9%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

print(bayesfactor_inclusion(bf))

Inclusion Bayes Factors (Model Averaged)

         P(prior) P(posterior) Inclusion BF
pert         0.60         1.00       322.52
cue          0.60         0.32        0.319
cue:pert     0.20         0.09        0.386

* Compared among: all models
*    Priors odds: uniform-equal

# Note: default multiplicity correction for pairs summary is Tukey HSD
print("----------------------------- EMMs -----------------------------")

[1] "----------------------------- EMMs -----------------------------"

emms_obj <- emmeans(temp_ANOVA, "pert")

NOTE: Results may be misleading due to involvement in interactions

print(pwpm(emms_obj, diffs = FALSE))

        curved   rot30
curved [0.451]  <.0001
rot30          [0.203]

Row and column labels: pert
Upper triangle: P values 
Diagonal: [Estimates] (emmean) 

# Effect size is Cohen's D
eff_size(emms_obj, edf = temp_ANOVA$anova_table$`den Df`[1], sigma = sigma(temp_ANOVA$lm))

 contrast       effect.size    SE df lower.CL upper.CL
 curved - rot30       0.956 0.235 77    0.487     1.42

Results are averaged over the levels of: cue 
sigma used for effect sizes: 0.2602 
Confidence level used: 0.95 

print("----------------------------- Full Tukey HSD -----------------------------")

[1] "----------------------------- Full Tukey HSD -----------------------------"

print(TukeyHSD(temp_ANOVA$aov))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = tmp_formula, data = dat.ret, contrasts = contrasts)

$pert
                   diff        lwr        upr    p adj
rot30-curved -0.2483742 -0.3635136 -0.1332348 5.03e-05

$cue
                       diff         lwr       upr     p adj
uncued-cued_tilt 0.05932141 -0.05581798 0.1744608 0.3081406

$`pert:cue`
                                        diff         lwr          upr     p adj
rot30:cued_tilt-curved:cued_tilt -0.20908884 -0.42515016  0.006972473 0.0614883
curved:uncued-curved:cued_tilt    0.09933909 -0.11672222  0.315400403 0.6240862
rot30:uncued-curved:cued_tilt    -0.18881487 -0.40228853  0.024658785 0.1017557
curved:uncued-rot30:cued_tilt     0.30842793  0.09236662  0.524489245 0.0019055
rot30:uncued-rot30:cued_tilt      0.02027397 -0.19319969  0.233747627 0.9945085
rot30:uncued-curved:uncued       -0.28815396 -0.50162762 -0.074680304 0.0036829

print(eff_size(emmeans(temp_ANOVA$aov, c("pert", "cue")),
  edf = temp_ANOVA$anova_table$`den Df`[1],
  sigma = sigma(temp_ANOVA$lm)
))

 contrast                           effect.size    SE df lower.CL upper.CL
 curved cued_tilt - rot30 cued_tilt      0.8036 0.323 77   0.1609    1.446
 curved cued_tilt - curved uncued       -0.3818 0.318 77  -1.0145    0.251
 curved cued_tilt - rot30 uncued         0.7257 0.318 77   0.0928    1.359
 rot30 cued_tilt - curved uncued        -1.1854 0.330 77  -1.8432   -0.528
 rot30 cued_tilt - rot30 uncued         -0.0779 0.313 77  -0.7002    0.544
 curved uncued - rot30 uncued            1.1075 0.325 77   0.4605    1.755

sigma used for effect sizes: 0.2602 
Confidence level used: 0.95 

1.6 Washout - Fit High Points (Starts)

p <- data_per_group %>%
  ggplot(
    aes(x = experiment, y = mean_high, colour = experiment)
  ) +
  theme_classic() +
  labs(
    x = per_group_x_labs,
    y = "Start Point of Decay"
  ) +
  # for the colour legend, only show the first 7
  # Note this doesn't work for the plotly plot
  guides(
    colour = guide_legend(override.aes = list(alpha = 1))
  ) +
  # set colour palette
  scale_colour_manual(values = pallete_list) +
  scale_fill_manual(values = pallete_list) +
  # add horizontal lines
  geom_hline(
    yintercept = c(0, 1), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "solid"
  ) +
  geom_hline(
    yintercept = c(0.5, 1.5), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "dashed"
  ) +
  # y-axis scale
  scale_y_continuous(
    limits = c(-0.5, 1.6),
    breaks = c(-0.5, 0, 0.5, 1, 1.5),
    labels = c(-0.5, 0, 0.5, 1, 1.5)
  ) +
  scale_x_discrete(
    labels = c(0, 0, 0, 0)
  ) +
  # remove all x axis labels
  theme(
    legend.position = "none",
    text = element_text(size = text_size)
  ) +
  # add data points
  geom_beeswarm(
    data = data_per_ppt,
    aes(
      y = exp_fit_N0
    ),
    alpha = 0.1,
    size = 1
  ) +
  geom_linerange(aes(
    ymin = mean_high - ci_high,
    ymax = mean_high + ci_high
  ), alpha = 0.5, lwd = 2) +
  geom_point()

# save
if (save_plots) {
  ggsave(
    p,
    filename = "../data/figs/paper_figs/fig3_washout_starts.pdf", device = "pdf",
    height = 2, width = per_group_fig_width
  )
}

p

1.6.1 ANOVA - 3 levels

(temp_ANOVA <- aov_car(
  exp_fit_N0 ~ pert * cue + Error(ppid),
  data_per_ppt,
  include_aov = TRUE
))

Contrasts set to contr.sum for the following variables: pert, cue

Anova Table (Type 3 tests)

Response: exp_fit_N0
    Effect    df  MSE         F  ges p.value
1     pert 1, 77 0.11 16.20 *** .174   <.001
2      cue 1, 77 0.11 18.31 *** .192   <.001
3 pert:cue 1, 77 0.11      0.83 .011    .364
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

print("Bayes anovas and bayes factors:")

[1] "Bayes anovas and bayes factors:"

print("")

[1] ""

bf <- anovaBF(exp_fit_N0 ~ pert * cue, data = data.frame(data_per_ppt), progress = FALSE)
print(bf)

Bayes factor analysis
--------------
[1] pert                  : 67.83776 ±0%
[2] cue                   : 151.9825 ±0%
[3] pert + cue            : 27387.72 ±2.26%
[4] pert + cue + pert:cue : 11973.49 ±0.94%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

print(bayesfactor_inclusion(bf))

Inclusion Bayes Factors (Model Averaged)

         P(prior) P(posterior) Inclusion BF
pert         0.60         1.00       171.82
cue          0.60         1.00       382.67
cue:pert     0.20         0.30         1.73

* Compared among: all models
*    Priors odds: uniform-equal

# Note: default multiplicity correction for pairs summary is Tukey HSD
print("----------------------------- EMMs -----------------------------")

[1] "----------------------------- EMMs -----------------------------"

emms_obj <- emmeans(temp_ANOVA, "pert")

NOTE: Results may be misleading due to involvement in interactions

print(pwpm(emms_obj, diffs = FALSE))

        curved   rot30
curved [0.527]  0.0001
rot30          [0.825]

Row and column labels: pert
Upper triangle: P values 
Diagonal: [Estimates] (emmean) 

# Effect size is Cohen's D
eff_size(emms_obj, edf = temp_ANOVA$anova_table$`den Df`[1], sigma = sigma(temp_ANOVA$lm))

 contrast       effect.size    SE df lower.CL upper.CL
 curved - rot30      -0.895 0.234 77    -1.36   -0.429

Results are averaged over the levels of: cue 
sigma used for effect sizes: 0.333 
Confidence level used: 0.95 

print("----------------------------- Full Tukey HSD -----------------------------")

[1] "----------------------------- Full Tukey HSD -----------------------------"

print(TukeyHSD(temp_ANOVA$aov))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = tmp_formula, data = dat.ret, contrasts = contrasts)

$pert
                  diff       lwr       upr     p adj
rot30-curved 0.3026107 0.1552576 0.4499637 0.0001054

$cue
                      diff       lwr       upr    p adj
uncued-cued_tilt 0.3174539 0.1701009 0.4648069 5.14e-05

$`pert:cue`
                                       diff         lwr       upr     p adj
rot30:cued_tilt-curved:cued_tilt 0.23030409 -0.04620673 0.5068149 0.1359571
curved:uncued-curved:cued_tilt   0.24906650 -0.02744432 0.5255773 0.0925397
rot30:uncued-curved:cued_tilt    0.61461041  0.34141122 0.8878096 0.0000005
curved:uncued-rot30:cued_tilt    0.01876241 -0.25774841 0.2952732 0.9979734
rot30:uncued-rot30:cued_tilt     0.38430632  0.11110712 0.6575055 0.0022792
rot30:uncued-curved:uncued       0.36554391  0.09234471 0.6387431 0.0040622

print(eff_size(emmeans(temp_ANOVA$aov, c("pert", "cue")),
  edf = temp_ANOVA$anova_table$`den Df`[1],
  sigma = sigma(temp_ANOVA$lm)
))

 contrast                           effect.size    SE df lower.CL upper.CL
 curved cued_tilt - rot30 cued_tilt     -0.6917 0.321 77   -1.331  -0.0523
 curved cued_tilt - curved uncued       -0.7480 0.322 77   -1.389  -0.1070
 curved cued_tilt - rot30 uncued        -1.8458 0.346 77   -2.535  -1.1568
 rot30 cued_tilt - curved uncued        -0.0563 0.316 77   -0.686   0.5734
 rot30 cued_tilt - rot30 uncued         -1.1542 0.326 77   -1.803  -0.5050
 curved uncued - rot30 uncued           -1.0978 0.325 77   -1.744  -0.4512

sigma used for effect sizes: 0.333 
Confidence level used: 0.95 

2 Alternate Washout Phases

These include the last trial of the Training phase for all Cued conditions.

2.1 Normalized Washout Phases

# rest of the exps
data_per_group <- omnibus_df %>%
  filter(
    experiment %in% c(
      "rot30_cued_tilt", "rot30_uncued",
      "curved_cued_tilt", "curved_uncued"
    ),
    trial_num >= 214,
    alt_washout_block == TRUE,
  ) %>%
  arrange(experiment, ppid, trial_num) %>%
  group_by(experiment, trial_num) %>%
  summarise(
    mean_deviation = mean(norm_throw_deviation),
    ci_deviation = vector_confint(norm_throw_deviation),
    .groups = "drop"
  )

# NOTE: this can't be combined with above since we are using nrow
data_per_group <- data_per_group %>%
  mutate(dummy_x = rep(1:40,
    length.out = nrow(data_per_group)
  ))

# add a dummy column with repeating sequence


# set up plot
p <- data_per_group %>%
  ggplot(
    aes(
      x = dummy_x, y = mean_deviation, colour = experiment
    )
  ) +
  theme_classic() +
  theme(
    legend.position = "none",
    text = element_text(size = text_size)
  ) +
  labs(
    x = "Trial",
    y = "Normalized Throw Angle"
  ) +
  # add horizontal lines
  geom_hline(
    yintercept = c(0, 1), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "solid"
  ) +
  geom_hline(
    yintercept = c(0.5, 1.5), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "dashed"
  ) +
  # set colour palette
  scale_colour_manual(values = pallete_list) +
  scale_fill_manual(values = pallete_list) + # scale axes
  scale_y_continuous(
    breaks = c(0, 0.5, 1, 1.5),
    labels = c(0, 0.5, 1, 1.5)
  ) +
  scale_x_continuous(
    breaks = curve_x_breaks,
    labels = curve_x_breaks
  )

# add confidence intervals and data points
for (uniq_exp in unique(data_per_group$experiment)) {
  # get the data for this block
  to_plot_data <- filter(data_per_group, experiment == uniq_exp)

  p <- p + geom_ribbon(
    data = to_plot_data,
    aes(
      ymin = mean_deviation - ci_deviation,
      ymax = mean_deviation + ci_deviation,
      fill = experiment
    ), colour = NA, alpha = 0.3
  ) + geom_line(
    data = to_plot_data
  )
}

# save
if (save_plots) {
  ggsave(
    p,
    filename = "../data/figs/paper_figs/fig3_curved_alt_washout.pdf", device = "pdf",
    height = 2, width = 2.5
  )
}

p

# ggplotly(p)

fitted_data <- fitted_alt_washout_df %>%
  filter(
    experiment %in% c(
      "rot30_cued_tilt", "rot30_uncued",
      "curved_cued_tilt", "curved_uncued"
    ),
    phase == "washout"
  ) %>%
  mutate(
    pert = factor(str_split_fixed(experiment, "_", 2)[, 1]),
    cue = factor(str_split_fixed(experiment, "_", 2)[, 2])
  )

p <- fitted_data %>%
  ggplot(
    aes(x = x + 1, y = y, colour = experiment)
  ) +
  theme_classic() +
  theme(
    legend.position = "none",
    text = element_text(size = text_size)
  ) +
  labs(
    x = "Trial",
    y = "Normalized Throw Angle"
  ) +
  # add horizontal lines
  geom_hline(
    yintercept = c(0, 1), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "solid"
  ) +
  geom_hline(
    yintercept = c(0.5, 1.5), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "dashed"
  ) +
  # set colour palette
  scale_colour_manual(values = pallete_list) +
  scale_fill_manual(values = pallete_list) + # scale axes
  scale_y_continuous(
    breaks = c(0, 0.5, 1, 1.5),
    labels = c(0, 0.5, 1, 1.5),
    limits = c(-0.25, 1.5)
  ) +
  scale_x_continuous(
    breaks = curve_x_breaks,
    labels = curve_x_breaks
  ) +
  geom_line()

# save
if (save_plots) {
  ggsave(
    p,
    filename = "../data/figs/paper_figs/fig3_curved_alt_washout_fitted_only.pdf", device = "pdf",
    height = 2, width = curve_fig_width
  )
}

p

2.2 Washout - Fit Learning Rates

data_per_ppt <- alt_washout_rates %>%
  filter(
    experiment %in% c(
      # "accel_cued_tilt", "accel_uncued",
      "rot30_cued_tilt", "rot30_uncued",
      "curved_cued_tilt", "curved_uncued"
    )
  ) %>%
  mutate(
    pert = factor(str_split_fixed(experiment, "_", 2)[, 1]),
    cue = factor(str_split_fixed(experiment, "_", 2)[, 2]),
    ppid = factor(ppid)
  )

data_per_group <- data_per_ppt %>%
  group_by(experiment) %>%
  summarise(
    mean_learning_rate = mean(exp_fit_lambda),
    ci_learning_rate = vector_confint(exp_fit_lambda),
    mean_high = mean(exp_fit_N0),
    ci_high = vector_confint(exp_fit_N0),
    .groups = "drop"
  )

p <- data_per_group %>%
  ggplot(
    aes(x = experiment, y = mean_learning_rate, colour = experiment)
  ) +
  theme_classic() +
  labs(
    x = per_group_x_labs,
    y = "Decay Rate"
  ) +
  # remove all x axis labels
  theme(
    legend.position = "none",
    text = element_text(size = text_size)
  ) +
  # colour legend settings
  guides(
    colour = guide_legend(override.aes = list(alpha = 1))
  ) +
  # set colour palette
  scale_colour_manual(values = pallete_list) +
  scale_fill_manual(
    values = pallete_list
  ) +
  # add horizontal lines
  geom_hline(
    yintercept = c(0, 1), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "solid"
  ) +
  geom_hline(
    yintercept = c(0.25, 0.5, 0.75), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "dashed"
  ) +
  # y-axis scale
  scale_y_continuous(
    # limits = c(-35, 5),
    breaks = c(0, 0.25, 0.5, 0.75, 1),
    labels = c(0, 0.25, 0.5, 0.75, 1)
  ) +
  scale_x_discrete(
    labels = c(0, 0, 0, 0)
  ) +
  # add data points
  geom_beeswarm(
    data = data_per_ppt,
    aes(
      y = exp_fit_lambda
    ),
    alpha = 0.1,
    size = 1
  ) +
  geom_linerange(aes(
    ymin = mean_learning_rate - ci_learning_rate,
    ymax = mean_learning_rate + ci_learning_rate
  ), alpha = 0.5, lwd = 2) +
  geom_point()

# save
if (save_plots) {
  ggsave(
    p,
    filename = "../data/figs/paper_figs/fig3_curved_alt_washout_learning_rates.pdf", device = "pdf",
    height = 2, width = per_group_fig_width
  )
}

p

(temp_ANOVA <- aov_car(
  exp_fit_lambda ~ pert * cue + Error(ppid),
  data_per_ppt,
  include_aov = TRUE
))

Contrasts set to contr.sum for the following variables: pert, cue

Anova Table (Type 3 tests)

Response: exp_fit_lambda
    Effect    df  MSE         F  ges p.value
1     pert 1, 77 0.05 44.94 *** .369   <.001
2      cue 1, 77 0.05      1.02 .013    .316
3 pert:cue 1, 77 0.05      0.58 .008    .447
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

print("Bayes anovas and bayes factors:")

[1] "Bayes anovas and bayes factors:"

print("")

[1] ""

bf <- anovaBF(exp_fit_lambda ~ pert * cue, data = data.frame(data_per_ppt), progress = FALSE)
print(bf)

Bayes factor analysis
--------------
[1] pert                  : 3524431   ±0%
[2] cue                   : 0.3210555 ±0.02%
[3] pert + cue            : 1199871   ±1.22%
[4] pert + cue + pert:cue : 487221.8  ±1.3%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

print(bayesfactor_inclusion(bf))

Inclusion Bayes Factors (Model Averaged)

         P(prior) P(posterior) Inclusion BF
pert         0.60         1.00     2.63e+06
cue          0.60         0.32        0.319
cue:pert     0.20         0.09        0.413

* Compared among: all models
*    Priors odds: uniform-equal

Pairwise comparisons:

# posthoc tests
# (m1 <- emmeans(temp_ANOVA, ~pert))
# print("Pairwise comparisons:")
# print(summary(pairs(m1), adjust = "tukey")) # Note: default multiplicity correction for pairs summary is Tukey HSD

# Note: default multiplicity correction for pairs summary is Tukey HSD
print("-------------------------- EMMs cue --------------------------")

[1] "-------------------------- EMMs cue --------------------------"

emms_obj <- emmeans(temp_ANOVA, "cue")

NOTE: Results may be misleading due to involvement in interactions

print(pwpm(emms_obj, diffs = FALSE))

          cued_tilt  uncued
cued_tilt   [0.406]  0.3159
uncued              [0.357]

Row and column labels: cue
Upper triangle: P values 
Diagonal: [Estimates] (emmean) 

# Effect size is Cohen's D
eff_size(emms_obj, edf = temp_ANOVA$anova_table$`den Df`[1], sigma = sigma(temp_ANOVA$lm))

 contrast           effect.size    SE df lower.CL upper.CL
 cued_tilt - uncued       0.224 0.223 77    -0.22    0.668

Results are averaged over the levels of: pert 
sigma used for effect sizes: 0.2183 
Confidence level used: 0.95 

print("-------------------------- EMMs perts --------------------------")

[1] "-------------------------- EMMs perts --------------------------"

emms_obj <- emmeans(temp_ANOVA, "pert")

NOTE: Results may be misleading due to involvement in interactions

print(pwpm(emms_obj, diffs = FALSE))

        curved   rot30
curved [0.544]  <.0001
rot30          [0.219]

Row and column labels: pert
Upper triangle: P values 
Diagonal: [Estimates] (emmean) 

# Effect size is Cohen's D
eff_size(emms_obj, edf = temp_ANOVA$anova_table$`den Df`[1], sigma = sigma(temp_ANOVA$lm))

 contrast       effect.size    SE df lower.CL upper.CL
 curved - rot30        1.49 0.253 77    0.987     1.99

Results are averaged over the levels of: cue 
sigma used for effect sizes: 0.2183 
Confidence level used: 0.95 

print("----------------------------- Full Tukey HSD -----------------------------")

[1] "----------------------------- Full Tukey HSD -----------------------------"

print(TukeyHSD(temp_ANOVA$aov))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = tmp_formula, data = dat.ret, contrasts = contrasts)

$pert
                   diff        lwr        upr p adj
rot30-curved -0.3254002 -0.4219955 -0.2288049     0

$cue
                        diff        lwr        upr     p adj
uncued-cued_tilt -0.04852349 -0.1451188 0.04807182 0.3203069

$`pert:cue`
                                        diff         lwr         upr     p adj
rot30:cued_tilt-curved:cued_tilt -0.36235677 -0.54361976 -0.18109378 0.0000077
curved:uncued-curved:cued_tilt   -0.08607912 -0.26734211  0.09518387 0.5991191
rot30:uncued-curved:cued_tilt    -0.37423308 -0.55332518 -0.19514099 0.0000030
curved:uncued-rot30:cued_tilt     0.27627765  0.09501466  0.45754064 0.0008105
rot30:uncued-rot30:cued_tilt     -0.01187631 -0.19096841  0.16721578 0.9981073
rot30:uncued-curved:uncued       -0.28815396 -0.46724606 -0.10906187 0.0003720

print(eff_size(emmeans(temp_ANOVA$aov, c("pert", "cue")),
  edf = temp_ANOVA$anova_table$`den Df`[1],
  sigma = sigma(temp_ANOVA$lm)
))

 contrast                           effect.size    SE df lower.CL upper.CL
 curved cued_tilt - rot30 cued_tilt      1.6601 0.343 77    0.976    2.344
 curved cued_tilt - curved uncued        0.3944 0.318 77   -0.239    1.027
 curved cued_tilt - rot30 uncued         1.7145 0.342 77    1.034    2.395
 rot30 cued_tilt - curved uncued        -1.2657 0.332 77   -1.927   -0.604
 rot30 cued_tilt - rot30 uncued          0.0544 0.312 77   -0.568    0.677
 curved uncued - rot30 uncued            1.3201 0.330 77    0.663    1.977

sigma used for effect sizes: 0.2183 
Confidence level used: 0.95 

2.3 Washout - Fit High Points (Starts)

p <- data_per_group %>%
  ggplot(
    aes(x = experiment, y = mean_high, colour = experiment)
  ) +
  theme_classic() +
  labs(
    x = per_group_x_labs,
    y = "Start Point of Decay"
  ) +
  # for the colour legend, only show the first 7
  # Note this doesn't work for the plotly plot
  guides(
    colour = guide_legend(override.aes = list(alpha = 1))
  ) +
  # set colour palette
  scale_colour_manual(values = pallete_list) +
  scale_fill_manual(values = pallete_list) +
  # add horizontal lines
  geom_hline(
    yintercept = c(0, 1), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "solid"
  ) +
  geom_hline(
    yintercept = c(0.5, 1.5), linewidth = 0.4,
    colour = "#CCCCCC", linetype = "dashed"
  ) +
  # y-axis scale
  scale_y_continuous(
    limits = c(-0.5, 1.6),
    breaks = c(-0.5, 0, 0.5, 1, 1.5),
    labels = c(-0.5, 0, 0.5, 1, 1.5)
  ) +
  scale_x_discrete(
    labels = c(0, 0, 0, 0)
  ) +
  # remove all x axis labels
  theme(
    legend.position = "none",
    text = element_text(size = text_size)
  ) +
  # add data points
  geom_beeswarm(
    data = data_per_ppt,
    aes(
      y = exp_fit_N0
    ),
    alpha = 0.1,
    size = 1
  ) +
  geom_linerange(aes(
    ymin = mean_high - ci_high,
    ymax = mean_high + ci_high
  ), alpha = 0.5, lwd = 2) +
  geom_point()

# save
if (save_plots) {
  ggsave(
    p,
    filename = "../data/figs/paper_figs/fig3_curved_alt_washout_starts.pdf", device = "pdf",
    height = 2, width = per_group_fig_width
  )
}

p

(temp_ANOVA <- aov_car(
  exp_fit_N0 ~ pert * cue + Error(ppid),
  data_per_ppt,
  include_aov = TRUE
))

Contrasts set to contr.sum for the following variables: pert, cue

Anova Table (Type 3 tests)

Response: exp_fit_N0
    Effect    df  MSE      F  ges p.value
1     pert 1, 77 0.14 3.73 + .046    .057
2      cue 1, 77 0.14   2.58 .032    .112
3 pert:cue 1, 77 0.14 6.10 * .073    .016
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

print("Bayes anovas and bayes factors:")

[1] "Bayes anovas and bayes factors:"

print("")

[1] ""

bf <- anovaBF(exp_fit_N0 ~ pert * cue, data = data.frame(data_per_ppt), progress = FALSE)
print(bf)

Bayes factor analysis
--------------
[1] pert                  : 1.034682  ±0.01%
[2] cue                   : 0.5958974 ±0.01%
[3] pert + cue            : 0.6489424 ±1.85%
[4] pert + cue + pert:cue : 2.308289  ±1%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

print(bayesfactor_inclusion(bf))

Inclusion Bayes Factors (Model Averaged)

         P(prior) P(posterior) Inclusion BF
pert         0.60         0.71         1.67
cue          0.60         0.64         1.16
cue:pert     0.20         0.41         2.82

* Compared among: all models
*    Priors odds: uniform-equal

Pairwise comparisons:

# posthoc tests
# (m1 <- emmeans(temp_ANOVA, ~pert))
# print("Pairwise comparisons:")
# print(summary(pairs(m1), adjust = "tukey")) # Note: default multiplicity correction for pairs summary is Tukey HSD

# Note: default multiplicity correction for pairs summary is Tukey HSD
print("-------------------------- EMMs cue --------------------------")

[1] "-------------------------- EMMs cue --------------------------"

emms_obj <- emmeans(temp_ANOVA, "cue")

NOTE: Results may be misleading due to involvement in interactions

print(pwpm(emms_obj, diffs = FALSE))

          cued_tilt  uncued
cued_tilt   [0.967]  0.1125
uncued              [0.834]

Row and column labels: cue
Upper triangle: P values 
Diagonal: [Estimates] (emmean) 

# Effect size is Cohen's D
eff_size(emms_obj, edf = temp_ANOVA$anova_table$`den Df`[1], sigma = sigma(temp_ANOVA$lm))

 contrast           effect.size    SE df lower.CL upper.CL
 cued_tilt - uncued       0.357 0.224 77  -0.0895    0.803

Results are averaged over the levels of: pert 
sigma used for effect sizes: 0.3737 
Confidence level used: 0.95 

print("-------------------------- EMMs perts --------------------------")

[1] "-------------------------- EMMs perts --------------------------"

emms_obj <- emmeans(temp_ANOVA, "pert")

NOTE: Results may be misleading due to involvement in interactions

print(pwpm(emms_obj, diffs = FALSE))

        curved   rot30
curved [0.821]  0.0572
rot30          [0.981]

Row and column labels: pert
Upper triangle: P values 
Diagonal: [Estimates] (emmean) 

# Effect size is Cohen's D
eff_size(emms_obj, edf = temp_ANOVA$anova_table$`den Df`[1], sigma = sigma(temp_ANOVA$lm))

 contrast       effect.size    SE df lower.CL upper.CL
 curved - rot30      -0.429 0.225 77   -0.877   0.0188

Results are averaged over the levels of: cue 
sigma used for effect sizes: 0.3737 
Confidence level used: 0.95 

print("----------------------------- Full Tukey HSD -----------------------------")

[1] "----------------------------- Full Tukey HSD -----------------------------"

print(TukeyHSD(temp_ANOVA$aov))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = tmp_formula, data = dat.ret, contrasts = contrasts)

$pert
                  diff          lwr       upr     p adj
rot30-curved 0.1612721 -0.004120548 0.3266647 0.0558383

$cue
                       diff        lwr        upr     p adj
uncued-cued_tilt -0.1308731 -0.2962657 0.03451958 0.1192056

$`pert:cue`
                                        diff         lwr         upr     p adj
rot30:cued_tilt-curved:cued_tilt -0.04475055 -0.35511304  0.26561194 0.9813751
curved:uncued-curved:cued_tilt   -0.33851140 -0.64887390 -0.02814891 0.0270589
rot30:uncued-curved:cued_tilt     0.02703251 -0.27961294  0.33367795 0.9955934
curved:uncued-rot30:cued_tilt    -0.29376085 -0.60412335  0.01660164 0.0702030
rot30:uncued-rot30:cued_tilt      0.07178306 -0.23486239  0.37842850 0.9270944
rot30:uncued-curved:uncued        0.36554391  0.05889847  0.67218935 0.0129087

print(eff_size(emmeans(temp_ANOVA$aov, c("pert", "cue")),
  edf = temp_ANOVA$anova_table$`den Df`[1],
  sigma = sigma(temp_ANOVA$lm)
))

 contrast                           effect.size    SE df lower.CL upper.CL
 curved cued_tilt - rot30 cued_tilt      0.1197 0.316 77   -0.510    0.750
 curved cued_tilt - curved uncued        0.9057 0.325 77    0.259    1.552
 curved cued_tilt - rot30 uncued        -0.0723 0.312 77   -0.695    0.550
 rot30 cued_tilt - curved uncued         0.7860 0.323 77    0.144    1.428
 rot30 cued_tilt - rot30 uncued         -0.1921 0.313 77   -0.815    0.431
 curved uncued - rot30 uncued           -0.9781 0.322 77   -1.620   -0.336

sigma used for effect sizes: 0.3737 
Confidence level used: 0.95